class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 3: Percepción remota, análisis masivo y GEE ###Procesamiento masivo de datos: Paralelización de rásters José A. Lastra<br> <a href="http://github.com/JoseLastra"> Github: JoseLastra</a><br> <a href="mailto:jose.lastra@pucv.cl"> jose.lastra@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Noviembre, 2023</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ - Paralelización: * LibrerÃa parallel * LibrerÃa foreach y doParallel - LibrerÃa terra: * app() * sapp() * lapply() ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> © Allison Horst ] --- ## Recordando algunos conceptos -- - Problemas de cálculo ([Jones, 2017](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html)): * **CPU**: proceso toma mucho tiempo de CPU. * **Memoria**: Demanda demasiada memoria. * **I/O**: toma mucho tiempo la lectura/escritura desde el disco. * **Network**: La red toma demasiado tiempo en transferir. -- - La paralelización viene a "solucionar" los problemas asociados a la CPU, dada la presencia de múltiples núcleos en los pc's modernos. --- ## Paralelización: Conceptos <center><img src="data:image/png;base64,#sock_fork.png" width="900px"/></center> <center>Fuente: Benito, 2021</center> --- ## foreach y doParallel: Ejemplo práctico con archivos ráster -- - Crearemos una lista de archivos **MOD13Q1** correspondientes al Ãndice NDVI desde el año 2000 hasta el 2021. -- - Ejecutaremos dos tareas simples: (1) cortar y (2) enmascarar datos a un área de estudio. -- - Usaremos los siguientes archivos: (1) **san_javier.gpkg**, (2) **MOD13Q1_dates_ts492.csv** y (3) **MOD13Q1_bands (carpeta)** ```r ## Listando archivos bandas <- list.files(path = 'MOD13Q1_bands/',pattern = '*.tif',full.names = T) ## archivo de fechas tabla <- read_csv('MOD13Q1_dates_ts492.csv') ## shape de San Javier shp <- vect('san_javier.gpkg') ``` --- ## foreach y doParallel: Ejemplo práctico con rásters -- - Preparación de clúster de procesamiento. ```r ## nombres de salida archivos nombres <- paste0('San_Javier_MOD13Q1/', tabla$dates,'_', tabla$ID,'_San_Javier.tif') ## Cluster set up nCluster <- detectCores(logical = F)-1 #Número de núcleos fÃsicos cl <- makeCluster(nCluster,type = 'PSOCK') #Configuración del cluster registerDoParallel(cl) #registrar para doParallel ``` --- ## foreach y doParallel: Ejemplo práctico con raster -- - Configuración de foreach y detención del clúster ```r foreach (i= 1:length(nombres), .packages = c('terra', 'magrittr')) %dopar% { r <- bandas[i] %>% rast() #crear banda en R m <- r %>% crop(shp) %>% mask(shp) # #cortar y enmascarar # guardar resultado writeRaster(m, filename = nombres[i], overwrite=T, datatype = 'INT2S') rm(r, m) #liberar memoria } stopCluster(cl) ``` -- - ¿Qué falló? --- ## foreach y doParallel: Ejemplo práctico con raster <center><img src="data:image/png;base64,#error_pointer.png" width="900px"/></center> -- - Esto lo podemos solucionar de varias maneras: - Utilizar **sf** para leer nuestro vector: `shp <- read_sf('archivo.gpkg')` - Pasar la ruta del archivo y leerlo en cada procesamiento: `ruta <- 'san_javier.gpkg'` y `shp <- vect(ruta)` dentro del **foreach** - Empaquetar nuestro SpatVector `shp <- vect('archivo.gpkg') %>% wrap()` y `vect(shp)` ó `unwrap(shp)` en nuestro **foreach** - O se puede recurrir a librerÃas como stars o similares que permitan serializar este tipo de procesos sin empaquetar los objetos. --- <center><img src="data:image/png;base64,#cpu.png" width="800px"/></center> <center>LabGRS, 2022</center> --- ## terra y la paralelización -- - Los ejemplos de paralelización vistos anteriormente, donde se envÃa información a cada nodo, no son del todo compatibles con los objetos de la librerÃa **terra** ([ver más, pág. 6](https://cran.r-project.org/web/packages/terra/terra.pdf)). -- - Podemos ajustar este comportamiento como ya vimos, sin embargo algunas cosas no pueden ser trabajadas de esta forma. -- - Por esta razón, terra provee funciones que permiten realizar operaciones empleando múltiples núcleos. Siendo **app()** la más común y simple de emplear. ```r ndvi_sj <- list.files(path = 'San_Javier_MOD13Q1/', pattern = '*.tif', full.names = T) %>% rast() ``` --- ## app() -- - Esta función nos permite operar en un **SpatRaster** multi o monobanda para aplicar una función normalmente de resumen (sum, min, max, median, etc.). También podemos operar en **SpatRasterDataset**. -- - **Importante**: dado la eficiencia computacional de **terra** muchas funciones resumen pueden ser aplicadas directamente con un funcionamiento óptimo. -- - Sin embargo, para archivos muy grandes o funciones complejas, puede ser útil la aplicación de **app()** ```r mean_noApp <- ndvi_sj %>% mean(na.rm = T) #simple mean_app <- app(ndvi_sj, fun = mean, na.rm = T, cores = 5) #5 cores ``` ``` ## Unit: seconds ## expr min lq mean median uq max neval ## simple_ver 1.244800 1.257324 1.312691 1.347865 1.353643 1.359821 5 ## cores_ver 1.228598 1.259428 1.284805 1.288690 1.312214 1.335096 5 ``` -- - **app()** opera los archivos de la misma forma que **apply()**: donde cada banda es una columna en la matriz y los pÃxeles son las filas. --- <img src="data:image/png;base64,#DIPGEOPR_03_8_files/figure-html/unnamed-chunk-7-1.png" width="100%" /> --- ## sapp() -- - Si lo que queremos es aplicar una función a cada banda dentro de nuestro **SpatRaster**, podemos emplear *sapp()*. ```r ndvi <- list.files(path = 'MOD13Q1_bands/', pattern = '*.tif', full.names = T) %>% rast() shp_v <- vect('san_javier.gpkg') #version simple ndvi_sj <- ndvi %>% crop(shp_v) %>% mask(shp_v) ``` ``` ## Unit: seconds ## expr min lq mean median uq max neval ## simple_process 2.000274 2.081524 2.174873 2.130812 2.21773 2.787479 30 ``` --- ## sapp() ```r # version con sapp() ## Creando la funcion crop_mask <- function(x, ...){ r <- x %>% crop(shp_v) %>% mask(shp_v) r } ##aplicacion ndvi_sj_sapp <- sapp(ndvi_sj, fun = crop_mask) ``` ``` ## Unit: seconds ## expr min lq mean median uq max neval ## sapp_process 5.247624 5.668062 5.855689 5.886459 5.989657 6.48664 5 ``` --- ## Observaciones -- - **sapp()** invoca a **sapply()** por detrás y aún no se ha implementado la paralelización para esta alternativa (**Información actualizada al 23 de Octubre de 2023**). -- - Tal cuál lo indica la documentación: *"(...) La paralelización es soportada, pero a menudo no ayuda, e incluso puede ser más lenta."* -- - Justificaciones para paralelizar con terra: * Disponemos de más de 8 núcleos para procesar. * Tenemos un archivo muy grande que trabajar. * O tenemos una función muy compleja (lenta) para aplicar sobre nuestros archivos. -- - Muchas veces el uso de **lapply()** puede ser suficiente. --- ## lapply() con terra -- - Tomemos nuestra lista de archivos originales y la función creada anteriormente para cortar y enmascarar. ```r ndvi_list <- list.files(path = 'MOD13Q1_bands/', pattern = '*.tif', full.names = T) #adaptando la función para lapply crop_mask <- function(x, ...){ x <- rast(x) #SpatRaster r <- x %>% crop(shp_v) %>% mask(shp_v) #aplicación de operación rm(x) # liberar memoria r #entregar objeto } ndvi_crop <- lapply(ndvi_list,FUN = crop_mask) ``` --- ## Recomendaciones -- - Siempre es mejor probar primero el código en un área pequeña antes de ejecutar el procesamiento de gran escala -- - Podemos usar un área de prueba transformada en un tibble e ir porbando pixel a pixel la función, facilitando la búsqueda de errores. -- - También podemos emplear otras estructuras de código como transformar nuestros objetos en listas y data frames para paralelizar empleando `future_lapply()` --- ## BibliografÃa complementaria - Bosco, J. (2018). "R Para Principantes". [Ver](https://bookdown.org/jboscomendoza/r-principiantes4/) - Cao, R. y Fernández, R. (2022). "Introducción al procesamiento en paralelo en R". En: "Técnicas de Remuestreo". [Ver](https://rubenfcasal.github.io/book_remuestreo/intro-hpc.html) - Gillespie, C., & Lovelace, R. (2021). "Efficient R programming: a practical guide to smarter programming." O'Reilly Media, Inc. [Ver](https://csgillespie.github.io/efficientR/) - Hijmans, R. J., Bivand, R., Forner, K., Ooms, J., Pebesma, E., & Sumner, M. D. (2022). Package ‘terra’. [Ver](https://brieger.esalq.usp.br/CRAN/web/packages/terra/terra.pdf) - Jones, M. (2017). "Quick Intro to Parallel Computing in R". [Ver](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html) - McCallum, E., & Weston, S. (2011). "Parallel R". O'Reilly Media, Inc. - Orellana, J. (2018). "HPC con R para investigadores". [Ver](https://bookdown.org/content/1498/) --- class: inverse middle 